Cell Rank - Trajectory Inference - Control Cohort¶

Goal: Mirror the trajectory analysis performed on the scCLEAN condition

1 - perform Cell Rank according to tutorials: https://cellrank.readthedocs.io/en/stable/index.html

2 - How many lineages are identified?

2 - Do those lineages correlate with coronary vs pulmonary artery?

In [54]:
import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import seaborn as sns
import scFates as scf

scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo", dpi_save=600)
cr.settings.verbosity = 2
In [95]:
#open the filtered batch corrected anndata matrix 
#samples underwent QC individually, and then aggregated using 'pegasus aggregate_matrix' CLI
#See notebook for identification of Force Directed Layout Trajectory graphs
adata = sc.read_h5ad('VSMC/PC_miQC/PC_miQC_subtypes_Control_VSMC_all_patients.h5ad')
adata
Out[95]:
AnnData object with n_obs × n_vars = 34890 × 16168
    obs: 'n_genes', 'n_counts', 'percent_mito', 'Channel', 'scale', 'louvain_labels', 'Tissue', 'doublet_score', 'pred_dbl', 'splice', 'value'
    var: 'featureid', 'n_cells', 'percent_cells', 'robust', 'highly_variable_features', 'mean', 'var', 'hvf_loess', 'hvf_rank'
    uns: 'PCs', '_attr2type', 'diffmap_evals', 'genome', 'louvain_labels_colors', 'louvain_resolution', 'modality', 'ncells', 'norm_count', 'pca', 'pca_features', 'stdzn_max_value', 'stdzn_mean', 'stdzn_std'
    obsm: 'X_diffmap', 'X_fle', 'X_pca', 'X_pca_harmony', 'X_phi', 'X_umap', 'diffmap_knn_distances', 'diffmap_knn_indices', 'pca_harmony_knn_distances', 'pca_harmony_knn_indices'
    varm: 'means', 'partial_sum'
    obsp: 'W_diffmap', 'W_pca_harmony'
In [33]:
adata.obs
Out[33]:
n_genes n_counts percent_mito Channel scale louvain_labels Tissue doublet_score pred_dbl splice value
barcodekey
P1_Coronary-AAACCCAAGTCTAGAA 4458 22340 7.488809 P1_Coronary 4.476276 1 Coronary 0.001757 False Yes 0
P1_Coronary-AAACCCACAGCACCCA 6294 40694 3.351845 P1_Coronary 2.457365 3 Coronary 0.001612 False Yes 0
P1_Coronary-AAACCCAGTTCCGCTT 3739 14339 3.528837 P1_Coronary 6.974473 1 Coronary 0.001371 False Yes 0
P1_Coronary-AAACGAAAGAAAGCGA 4753 26268 4.301812 P1_Coronary 3.807058 2 Coronary 0.002118 False Yes 0
P1_Coronary-AAACGAAAGGTGGGTT 3288 13315 3.604957 P1_Coronary 7.510327 3 Coronary 0.001097 False Yes 0
... ... ... ... ... ... ... ... ... ... ... ...
P4_Pulmonary-TTTGGTTAGACTTCGT 6644 61325 1.125153 P4_Pulmonary 1.630683 2 Pulmonary 0.001096 False Yes 0
P4_Pulmonary-TTTGTTGAGCATCAAA 6367 40162 2.445097 P4_Pulmonary 2.490040 2 Pulmonary 0.001804 False Yes 0
P4_Pulmonary-TTTGTTGAGGTCTGGA 5293 35347 3.010156 P4_Pulmonary 2.829254 1 Pulmonary 0.004670 False Yes 0
P4_Pulmonary-TTTGTTGGTGACCTGC 5182 30960 3.003876 P4_Pulmonary 3.229974 1 Pulmonary 0.000443 False Yes 0
P4_Pulmonary-TTTGTTGTCTACCTTA 7676 61309 2.089416 P4_Pulmonary 1.631188 2 Pulmonary 0.000931 False Yes 0

34890 rows × 11 columns

In [34]:
#count the proper number of tissues to ensure the number of colors matches
adata.obs['Tissue'].value_counts()
Out[34]:
Pulmonary    19297
Coronary     15593
Name: Tissue, dtype: int64
In [35]:
sc.tl.embedding_density(adata, basis='fle', groupby='Tissue')
In [36]:
sc.set_figure_params(figsize=(6,3),frameon=False,dpi_save=300)
sc.pl.embedding_density(
    adata,
    basis='fle',
    key='fle_density_Tissue',
    color_map='GnBu',
    save='Density_Control_coronary_pulmonary.png')
WARNING: saving figure to file figures/fle_density_Tissue_Density_Control_coronary_pulmonary.png
/Users/jbezney/opt/anaconda3/envs/CELL_RANK/lib/python3.9/site-packages/scanpy/_settings.py:447: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`
  IPython.display.set_matplotlib_formats(*ipython_format)
/Users/jbezney/opt/anaconda3/envs/CELL_RANK/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:444: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  pl.colorbar(cax, ax=ax, pad=0.01, fraction=0.08, aspect=30)
In [37]:
adata.uns['Tissue_colors']= np.array(['#BA0900', '#809693'], dtype=object)
In [38]:
#need to use the harmony batch integrated matrix
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
In [39]:
# see Cell Rank tutorial "CellRank beyond RNA velocity"
# we're copying our `.X` matrix into the layers because that's where `scv.tl.moments`
# expects to find counts for imputation
adata.layers["spliced"] = adata.X
adata.layers["unspliced"] = adata.X
#30 and 30 
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
computing neighbors
    finished (0:00:02) --> added 
    'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
    finished (0:00:57) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
In [40]:
from cellrank.tl.kernels import CytoTRACEKernel

ctk = CytoTRACEKernel(adata)
Computing CytoTRACE score with `16168` genes
Adding `adata.obs['ct_score']`
       `adata.obs['ct_pseudotime']`
       `adata.obs['ct_num_exp_genes']`
       `adata.var['ct_gene_corr']`
       `adata.var['ct_correlates']`
       `adata.uns['ct_params']`
    Finish (0:00:18)
In [73]:
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", color="Tissue", legend_loc="right")
In [74]:
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", c="louvain_labels", legend_loc="right")
In [75]:
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", c="ct_pseudotime", legend_loc="right",cmap='magma',
              title='10x-V3: CytoTRACE Pseudotime', dpi=300, save='VSMC/PC_Figures/Control_Pseudotime.png')
saving figure to file VSMC/PC_Figures/Control_Pseudotime.png
In [44]:
ctk.compute_transition_matrix(threshold_scheme="soft", nu=0.3)
Computing transition matrix based on `ct_pseudotime`
100%|███████████████████████████████████████████████████████████████████████████████████████████████| 34890/34890 [00:07<00:00, 4864.46cell/s]
    Finish (0:00:17)
Out[44]:
<CytoTRACEKernel>
In [45]:
from cellrank.tl.estimators import GPCCA

g_fwd = GPCCA(ctk)
print(g_fwd)
GPCCA[n=34890, kernel=<CytoTRACEKernel[dnorm=False, scheme=soft, b=10.0, nu=0.3]>]
In [72]:
g_fwd.compute_schur(n_components=25, alpha=0.2)
g_fwd.plot_spectrum(real_only=True, title='10x-V3: Eigengap after 2 Eigenvalues', dpi=300,show_all_xticks=False, 
                   save='VSMC/Figures/Schur_spectrum_control.png')
Computing Schur decomposition
Mat Object: 1 MPI processes
  type: seqdense
9.9999999999999978e-01 1.5104703446893388e-03 -1.0589381352350191e-03 2.7950273950795579e-02 5.7465896995287716e-03 -2.6981007403120198e-02 -1.3253115379381099e-02 -6.8190724781176497e-04 -2.6218359832608572e-02 -3.3066354637139618e-02 -5.6385359598725114e-03 1.2441219204278382e-02 1.7873711968797199e-02 -2.7108550577416082e-03 6.9743318436813308e-03 -2.4262715375852124e-02 -7.0706238090262205e-02 -2.5899948554389304e-02 -5.7455826297393542e-02 -9.5251526184787544e-04 6.2442456244687344e-03 -1.7232947388803023e-02 8.3140904450904451e-03 1.0946223742786766e-02 -2.1585266902328183e-02 1.2795137726302936e-02 2.4070688306005745e-02 -2.7892400598665958e-02 
0.0000000000000000e+00 9.9627336204564065e-01 5.4521970128935350e-03 9.8249454254520645e-04 -2.6367971432676867e-02 -1.0192294136377173e-02 2.4376187853096100e-04 4.5190890936599318e-03 -4.6776283811514763e-03 -1.9275232991329244e-02 6.5829359852476732e-03 -9.0155665327898236e-03 5.7279085374940559e-04 1.8412026023431481e-03 2.3991022827994819e-02 5.6427958936241017e-03 -3.9958310628242268e-03 -6.2546309197645993e-03 -2.5892796843553598e-02 -2.3725333701813622e-02 -2.1431039784099909e-02 -1.5122492873592554e-02 -1.7984925923894343e-04 3.8945076407209691e-03 -1.1261512805184245e-02 1.3967370317709480e-02 4.0139917876200592e-03 -1.1985099383889086e-02 
0.0000000000000000e+00 0.0000000000000000e+00 9.8324778875787489e-01 -8.7536686124265232e-04 -7.4549290234413364e-03 -1.8518199376798053e-03 5.8744103599715823e-04 1.9050220793504818e-03 9.2603528888584501e-04 -4.8631991096830130e-03 3.4140223575729383e-03 -6.5934815214324935e-03 -1.0810547788351638e-03 1.4496609039301213e-03 7.3787955764424342e-03 2.7261606675177451e-03 3.2964656703700426e-03 -1.1114521927562715e-04 -5.0941254438301520e-03 -1.1629490234482219e-02 -8.5878277955733981e-03 -4.4138513740445235e-03 -1.0870281142966222e-03 5.0709826635106597e-04 -4.7832693433645527e-03 5.2446725314915623e-03 5.0538103758852350e-04 -1.4454581112029679e-04 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.7281118371582853e-01 -3.4243483946466316e-03 5.6100879434757226e-03 1.1458710087541700e-03 9.9169271132793467e-03 1.4769863460627279e-02 -3.2024693464858776e-02 1.2453537831746632e-02 -7.8050864185200981e-03 -5.6355978550003596e-03 -3.1423887825281051e-03 3.3412439886181693e-02 -2.3104836168758941e-02 -1.4954065828873250e-02 8.0858369544636460e-03 -2.2799534812089885e-02 -1.8261462149155934e-03 -1.7462234997318536e-03 -1.5881001400214204e-02 2.2256300629418285e-02 1.9404651083168405e-02 -6.0907301755591080e-03 -1.2596475454979378e-02 5.0275210019492507e-03 2.7542019367107911e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.6172932137508682e-01 -8.7845109993197369e-03 -6.1183234209011613e-03 2.7683486888075357e-03 -4.8833301761978928e-04 -5.1794753792574078e-03 9.7072073606289423e-03 -5.6372498432818628e-02 -1.1814715755563980e-03 1.0468379782087318e-03 -3.8194964193303651e-03 7.9853912674039528e-03 9.1965814837766113e-03 -2.0766809648706148e-03 8.7569565838912910e-03 -1.0971524871770765e-02 -3.3215076391372389e-02 -3.8349354240162656e-02 2.6253722125541643e-02 4.3057621918894222e-02 3.0019044654752032e-02 9.7806195813078048e-03 1.1206195923855113e-02 -1.3101122027606577e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.4946879335557066e-01 -1.7090980012679654e-04 -1.0245330396447379e-02 -1.8195505975496695e-02 3.6766217943761594e-02 -1.0713035048589626e-03 -1.0706382856721553e-02 2.9890276086685909e-02 3.8056981527571379e-02 -2.2918869448018173e-02 2.1869217212117480e-02 6.1331634767788062e-03 -2.2821924301566973e-02 4.7976218655189085e-03 7.6767847554851387e-03 -2.0634331301406090e-02 2.1577285663641862e-02 1.7888410356856652e-02 -1.5174625671269159e-02 2.8892934827294404e-03 -8.0278714961815878e-03 3.0097383206820566e-02 2.2826542341472869e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.4855679846759411e-01 -1.2862596905909933e-02 -6.2428715908546069e-04 6.2736654924203670e-03 -6.9806715960521997e-03 -4.7159087161985176e-03 -3.9930468894956847e-03 -1.0789901096621040e-02 -5.2467957874303479e-02 -2.0395440728660508e-02 -5.4624990872841356e-03 -7.0795431961561957e-04 -9.6168405801234022e-03 -1.4203691296345612e-02 -7.9488779336637800e-03 4.9613427323571971e-03 4.4431671382889198e-03 1.9096326291047131e-02 -1.8294854954040707e-02 -3.9929091011835002e-03 -8.0514470496648588e-03 4.0451953589315505e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3259396545442486e-01 5.0664225718467680e-03 -5.3226975367169750e-03 3.1265498405836883e-03 3.0114429570133245e-03 -1.2228612209173356e-03 -5.7086896471761272e-03 3.7555447879614318e-02 -1.3640248583780384e-02 -3.7023152532333309e-03 7.1199075047312117e-03 -1.3515875430544576e-02 -8.3177962296825548e-03 2.0319152375604558e-02 6.1471860783078036e-03 6.7263133994036068e-03 3.5125755728424774e-03 -4.7053449032174420e-03 8.1749609431236245e-03 2.8391250885571758e-02 3.8228862708859911e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3158694969544364e-01 2.2242113804779071e-02 -1.7787193291740753e-02 -7.2900936945398259e-03 2.0269979633625746e-02 -6.6844659827845331e-03 -9.0457484384050930e-04 -8.5817263543636335e-03 1.3606310431806688e-02 -1.5827126767738290e-02 -2.8581881263804353e-02 1.6613578669786144e-02 -3.1633676743946244e-02 6.6038715541166498e-02 1.6471657918295527e-02 -1.0278054871803461e-02 1.8014980926325704e-02 -1.3560728951056034e-03 4.4823170021986429e-04 3.0886390633533153e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3091934441598767e-01 1.4223547926352253e-03 -7.4747519181522131e-04 -4.9652263254875553e-02 5.3607535341217386e-02 -2.2779534248062453e-02 -3.3482297464046587e-02 6.4111895820737980e-02 -1.0135346462070391e-02 -1.0861775008494346e-02 3.2400073128244139e-02 -2.2159238280316892e-02 1.1970995344664784e-02 -3.4941868871880916e-02 1.1929969989059696e-02 1.3990175813765771e-02 -4.7129354204493837e-02 9.1803260027871966e-03 5.2856007136041962e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2991898979791554e-01 1.3587882449117615e-03 1.1393210808298901e-02 -3.2747658221302448e-02 5.7516008468623736e-03 -1.8106653729116958e-02 -1.9787651274014297e-02 4.1712972130151911e-02 -9.4350303229051528e-03 -4.8614276959275384e-03 1.1261269689240464e-02 -4.4003551103625679e-03 -3.9008429620648891e-02 2.5707423059467934e-02 1.8671712238869488e-03 4.0434754961587237e-02 -3.2634577319055787e-02 -5.6334776181945850e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2434023999884818e-01 8.4716453730158178e-03 -1.3527031095716900e-02 1.1060045157607910e-02 7.4876374524998580e-03 -1.1551858947957062e-02 8.8025999974750500e-03 -4.2019534810376477e-03 -1.1374750256859918e-02 -2.2527009073594371e-02 -1.9150335689612559e-02 2.6264754049791554e-02 2.9646350339642732e-02 1.1523702629486302e-02 3.4045959764913611e-03 -1.5961648012890284e-02 -1.7314801715350974e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.1010320312826287e-01 -9.3367368129166940e-03 -1.3417063701558190e-03 5.3608466735433249e-02 -6.2745903247317311e-02 -4.2472737402125184e-02 -1.5292080723881928e-02 3.9842058609974070e-04 -2.5047939541090559e-02 -1.0474608058555660e-02 -2.6229788129411971e-02 4.6616788094295536e-03 -1.2815713736169653e-03 -6.0364812103400542e-02 -7.3459676396517961e-03 2.9940608194055512e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.0468483889695872e-01 -2.1542281197337955e-02 2.1786335835019113e-02 2.0392537366547787e-02 -1.5560562422447339e-02 2.4914337227839884e-03 3.7317987348866695e-03 1.8794826479532879e-02 -1.8739265637647724e-02 -5.1018114331322417e-02 -8.5335098852121691e-03 3.3761831483937361e-02 -1.3001025784296417e-03 -3.9055727948641213e-02 1.6685694701737919e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.0373563739105400e-01 -1.0288112854357262e-02 -3.1813185697879033e-02 -9.1366958887808009e-03 -4.0032634892655204e-02 -2.6667754453496686e-02 2.4844528849893847e-02 2.0599116213199879e-02 7.1591493556318183e-03 -8.2513569852320275e-03 -7.3159237525703311e-04 -4.2282238112723546e-02 1.3276309224545799e-02 2.6528015054847393e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8788208359093845e-01 1.3198667006794949e-03 2.3096728800521176e-03 5.6197966207582233e-03 -6.0187051411862624e-03 3.2586029140755805e-02 -3.2280725891437663e-02 -1.1619993121017169e-02 -2.7839459884860100e-03 -5.0938190149552309e-03 4.3674424707526881e-02 -1.3045673615349831e-02 4.5703692167642126e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8540604290006852e-01 -5.2134541440334373e-03 5.4269199247401348e-03 7.6339117295173503e-03 -2.6841612115697045e-03 1.8676373216803252e-02 -1.3716188320876102e-02 -7.7179079727982035e-03 -1.1871583928758622e-02 -1.1924843370880008e-02 2.1031912642256220e-02 -3.8951642414570440e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8281878410465775e-01 5.7374980637376477e-04 1.2291501419757667e-03 -4.5234218247302855e-03 1.3732616690500144e-02 -2.2669615996593024e-02 2.6285271741407368e-02 -1.8175093128818222e-03 1.2755980199534167e-02 1.6738455944954024e-02 1.0072490460355597e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7476711397366180e-01 -1.3366865117628298e-02 -1.8867374836053711e-03 1.9918392176175529e-02 3.5547592589591425e-03 -1.4017553556386293e-02 2.7515008565390599e-02 -2.2110532635193358e-02 -1.3644807946228297e-02 -9.7217976505553541e-04 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7333096535089350e-01 -1.0863173778571196e-02 -1.1433120341631534e-02 1.1776996810847274e-02 2.2805372551047851e-02 9.9717020028976290e-03 -3.1145766570008643e-03 1.0393064067801428e-02 1.4518733214855936e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7042956314455144e-01 -2.2648521206044311e-02 1.7114528820417699e-02 3.0723547412192902e-02 1.7062258250441002e-02 2.2879725100850734e-02 -1.0139789902297101e-02 -1.8119897238074605e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.6016086987201523e-01 1.7539035439345905e-02 7.1252777802771200e-03 6.0251560770741138e-03 -3.3921803409182480e-02 8.5303425146413935e-03 4.1715796998404650e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.5581726137817804e-01 -2.2127911254151956e-03 -2.1804594232151260e-02 -1.4663114190629576e-02 2.3199767463206001e-02 4.1731446761738199e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.5324121383954465e-01 -1.2541240200928335e-02 1.6663315099352170e-02 -1.2834388874028175e-02 -2.8735122943675678e-02 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4482251210635406e-01 3.2124717459147850e-03 -1.2006066397247104e-02 -2.7960591408977150e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4378110563514586e-01 3.2078737172057856e-03 -1.8088168604686499e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4294492729348747e-01 -2.4050511964898632e-04 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4188552489685931e-01 
Adding `adata.uns['eigendecomposition_fwd']`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:02)

According to the gap in the real component of the eigenvalues, two mcarostates should be computed

In [98]:
#cluster_key tries to associated the names of macrostates with cluster labels, 
#or in this case, tissue origin
g_fwd.compute_macrostates(n_states=2, cluster_key="Tissue")
g_fwd.plot_macrostates(
    discrete=True, legend_loc="right", size=100, basis="fle", title='scCLEAN Macrostates: 2',
    figsize=[8,8], dpi=400)
Computing `2` macrostates
Adding `.macrostates`
       `.macrostates_memberships`
       `.coarse_T`
       `.coarse_initial_distribution
       `.coarse_stationary_distribution`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:00)
In [99]:
#the course grained transition matrix is stable at 1
#that means there is a clear defined seperation of terminal states
g_fwd.plot_coarse_T(figsize=[5,6], title='10x-V3: Course-grained Transition Matrix', 
                  save='VSMC/PC_Figures/Course_Grained_Control.png', dpi=300)

the higher the stationary distance the more terminal end point in pseudotime

Stationary distance is the stationary distirbution of the transition matrix

The smaller the value the shorter the differentiation path

In [100]:
g_fwd.plot_macrostates(discrete=True, same_plot=True, basis="fle", size=50)
In [101]:
g_fwd.plot_macrostates(same_plot=False, basis='fle', 
                       title=['scCLEAN Macrostate 1','scCLEAN Macrostate 2'],
    figsize=[8,6], size=100, fontsize=18,  dpi=400)
In [102]:
#add a new observale into adata.obs defining the terminal states
g_fwd.set_terminal_states_from_macrostates(names=['Coronary',
                                                  'Pulmonary'])
Adding `adata.obs['terminal_states']`
       `adata.obs['terminal_states_probs']`
       `.terminal_states`
       `.terminal_states_probabilities`
       `.terminal_states_memberships
    Finish`
In [103]:
g_fwd.rename_terminal_states({'Pulmonary':'Lineage_1', 'Coronary':'Lineage_2'})
In [104]:
#this computes the fate probabilities toward those terminal states for every cell
g_fwd.compute_absorption_probabilities( time_to_absorption='all')
g_fwd.plot_absorption_probabilities(same_plot=False, size=50, basis="fle")
Computing absorption probabilities
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00,  3.12/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  1.49/s]
Adding `adata.obsm['to_terminal_states']`
       `adata.obsm['absorption_times_fwd']`
       `.absorption_probabilities`
       `.absorption_times`
    Finish (0:00:01)

Now it makes sense why the coronary macrostate had a stationary distance of 0

There are only a handful of cells differentiating toward this terminal state

99% of cells are transitioning toward macrostate associated with lineage 2

In [83]:
#check the probabilities of each cell moving to each fate 
lineage = pd.DataFrame(adata.obsm['to_terminal_states'])
lineage['Lineage_2'] = lineage[0]>0.5
lineage['Lineage_1'] = lineage[1]>0.5
lineage
Out[83]:
0 1 Lineage_2 Lineage_1
0 1.323540e-12 0.999660 False True
1 2.114814e-12 0.999660 False True
2 9.830771e-13 0.999660 False True
3 4.422072e-12 0.999660 False True
4 1.075595e-12 0.999660 False True
... ... ... ... ...
34885 8.360096e-10 0.999660 False True
34886 5.840236e-07 0.999660 False True
34887 1.398509e-09 0.999660 False True
34888 1.321043e-09 0.999660 False True
34889 1.363352e-06 0.999659 False True

34890 rows × 4 columns

In [84]:
#only 50 cells are differentiation toward lineage 1
lineage['Lineage_2'].value_counts()
Out[84]:
False    34840
True        50
Name: Lineage_2, dtype: int64
In [97]:
adata.obs['Tissue'].value_counts()
Out[97]:
Pulmonary    19297
Coronary     15593
Name: Tissue, dtype: int64
In [105]:
g_fwd.plot_absorption_probabilities(same_plot=False, size=50, 
                                    basis="fle", rescale_color=[0,1.001],
                                   title=['10x-V3 Lineage 2: Absorption Probability','10x-V3 Lineage 1: Absorption Probability'],
    figsize=[8,6], fontsize=18,  dpi=400, save='VSMC/PC_Figures/Control_absorption_probabilities.png')
saving figure to file VSMC/PC_Figures/Control_absorption_probabilities.png
In [85]:
#use any kwargs from https://scvelo.readthedocs.io/scvelo.pl.scatter/#scvelo.pl.scatter
g_fwd.plot_absorption_probabilities(same_plot=False, color='Tissue', mode='time', time_key='ct_pseudotime', 
                                    size=25, basis="fle", alpha=0.5, figsize=[8,8], fontsize=18,  dpi=400,
                                    xlabel='Pseudotime', ylabel='Probability',xlim=[-0.05,1.05],
                                   title=['10x-V3 Lineage 2: Absorption Probability','10x-V3 Lineage 1: Absorption Probability'],
                                   save='VSMC/PC_Figures/Control_probability_Coronary_Pulmonary.png')
saving figure to file VSMC/PC_Figures/Control_probability_Coronary_Pulmonary.png

Similar to scCLEAN condition, it is identifying that one lineage is capturing features specific to coronary

Except the amount of cells it can identify as probably differentiating toward that lineage is negligible

In [86]:
scv.pl.scatter(adata, basis="fle", c="terminal_states", size=100, legend_loc='on data',
               title='10x-V3 Terminal States', dpi=300, save='VSMC/PC_Figures/Control_1lineage_Terminal_States.png')
saving figure to file VSMC/PC_Figures/Control_1lineage_Terminal_States.png

AUC - Classification¶

What is the classification potential for these two lineage distinctions?

In [87]:
#read terminal state probabilities into dataframe to determine correlation
terminal_df = pd.DataFrame(adata.obsm['to_terminal_states'])
terminal_df = terminal_df.rename(columns={0: "Lineage 1", 1: "Lineage 2"})
terminal_df['barcode']=adata.obs.index.values
terminal_df = terminal_df.set_index('barcode')
terminal_df['Origin']=adata.obs['Tissue']
terminal_df
Out[87]:
Lineage 1 Lineage 2 Origin
barcode
P1_Coronary-AAACCCAAGTCTAGAA 1.323540e-12 0.999660 Coronary
P1_Coronary-AAACCCACAGCACCCA 2.114814e-12 0.999660 Coronary
P1_Coronary-AAACCCAGTTCCGCTT 9.830771e-13 0.999660 Coronary
P1_Coronary-AAACGAAAGAAAGCGA 4.422072e-12 0.999660 Coronary
P1_Coronary-AAACGAAAGGTGGGTT 1.075595e-12 0.999660 Coronary
... ... ... ...
P4_Pulmonary-TTTGGTTAGACTTCGT 8.360096e-10 0.999660 Pulmonary
P4_Pulmonary-TTTGTTGAGCATCAAA 5.840236e-07 0.999660 Pulmonary
P4_Pulmonary-TTTGTTGAGGTCTGGA 1.398509e-09 0.999660 Pulmonary
P4_Pulmonary-TTTGTTGGTGACCTGC 1.321043e-09 0.999660 Pulmonary
P4_Pulmonary-TTTGTTGTCTACCTTA 1.363352e-06 0.999659 Pulmonary

34890 rows × 3 columns

In [88]:
adata.obs['Pulmonary_lin_abs_prob']=terminal_df['Lineage 1']
adata.obs['Coronary_lin_abs_prob']=terminal_df['Lineage 2']
In [89]:
#change the identity state into a binary filter 
terminal_df['P_identity'] = 0
terminal_df.loc[terminal_df['Origin'] == 'Pulmonary', 'P_identity'] = 1

terminal_df['C_identity'] = 0
terminal_df.loc[terminal_df['Origin'] == 'Coronary', 'C_identity'] = 1

#create new dataframe to generate the heatmap
heat = pd.DataFrame()
heat['Lineage 1'] = terminal_df['Lineage 2']
heat['Pulmonary Identity'] = terminal_df['P_identity']
heat['Lineage 2'] = terminal_df['Lineage 1']
heat['Coronary Identity'] = terminal_df['C_identity']
heat = heat.reset_index(drop=True)
heat
Out[89]:
Lineage 1 Pulmonary Identity Lineage 2 Coronary Identity
0 0.999660 0 1.323540e-12 1
1 0.999660 0 2.114814e-12 1
2 0.999660 0 9.830771e-13 1
3 0.999660 0 4.422072e-12 1
4 0.999660 0 1.075595e-12 1
... ... ... ... ...
34885 0.999660 1 8.360096e-10 0
34886 0.999660 1 5.840236e-07 0
34887 0.999660 1 1.398509e-09 0
34888 0.999660 1 1.321043e-09 0
34889 0.999659 1 1.363352e-06 0

34890 rows × 4 columns

In [90]:
#plot the lineage probability compressed along the pseudotime x-axis
fig = plt.figure()

fig.set_figheight(8)
fig.set_figwidth(6)
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax = sc.pl.violin(adata, 'Coronary_lin_abs_prob', groupby='Tissue', stripplot=True, jitter=0.1, show=False, ax=ax)
ax.set_title('Lineage 2')
ax.set_ylabel('Absorption Probability')
ax.set_ylim(0, 1.05)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

ax2 = sc.pl.violin(adata, 'Pulmonary_lin_abs_prob', groupby='Tissue', stripplot=True, jitter=0.1, show=False, ax=ax2)
ax2.set_title('Lineage 1')
ax2.set_ylabel('Absorption Probability')
ax2.set_ylim(0, 1.05)

ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

fig.savefig('VSMC/PC_Figures/Control_Violin_Lineages.png', dpi=400)
In [91]:
lineage2 = heat[['Lineage 1','Pulmonary Identity']]

fig,ax = plt.subplots(figsize=(3,8))
ax = sns.heatmap(lineage2, cmap='Greys', yticklabels = 5000, 
                 )

plt.xticks(rotation = 30)
# set y-axis label
ax.set_ylabel("Cells",fontsize=12)
ax.figure.axes[-1].set_ylabel('Absorption Probability', size=12, rotation=270, labelpad=20)
fig.autofmt_xdate()

plt.title('Lineage 1 = Pulmonary Lineage', y=1.01, fontsize=12)

# save the plot as a file
fig.savefig('VSMC/PC_Figures/Control_Heatmap_Pulmonary.jpg', dpi=300, bbox_inches='tight')

plt.show()
In [92]:
lineage1 = heat[['Lineage 2','Coronary Identity']]

fig,ax = plt.subplots(figsize=(3,8))
ax = sns.heatmap(lineage1, cmap='Reds', yticklabels = 5000, 
                 )

plt.xticks(rotation = 30)
# set y-axis label
ax.set_ylabel("Cells",fontsize=12)
ax.figure.axes[-1].set_ylabel('Absorption Probability', size=12, rotation=270, labelpad=20)
fig.autofmt_xdate()

plt.title('Lineage 2 = Coronary Lineage', y=1.01, fontsize=12)

# save the plot as a file
fig.savefig('VSMC/PC_Figures/Control_Heatmap_Coronary.jpg', dpi=300, bbox_inches='tight')

plt.show()
In [93]:
#this is for the pulmonary lineage alone exclusively 
score = np.array(terminal_df['Lineage 2'])
y = np.array(terminal_df['P_identity'])

# false positive rate
fpr = []
# true positive rate
tpr = []
# Iterate thresholds from 0.0, 0.01, ... 1.0
thresholds = np.arange(0.0, 1.01, .01)

# get number of positive and negative examples in the dataset
P = sum(y)
N = len(y) - P

# iterate through all thresholds and determine fraction of true positives
# and false positives found at this threshold
for thresh in thresholds:
    FP=0
    TP=0
    for i in range(len(score)):
        if (score[i] > thresh):
            if y[i] == 1:
                TP = TP + 1
            if y[i] == 0:
                FP = FP + 1
    fpr.append(FP/float(N))
    tpr.append(TP/float(P))

import sklearn.metrics as metrics
roc_auc = metrics.auc(fpr, tpr)

fig, ax = plt.subplots()
fig.set_figheight(8)
fig.set_figwidth(6)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)


plt.title('ROC: Lineage 1 = Pulmonary Lineage')
plt.plot(fpr, tpr, 'gray', label = 'AUC = %0.2f' % roc_auc)
plt.plot([0, 1], [0, 1],'black',linestyle='dashed', label = 'Random Classifier')
plt.legend(loc = 'lower right')
plt.xlim([-0.02, 1])
plt.ylim([0, 1.05])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')

plt.savefig('VSMC/PC_Figures/Control_ROC_pulmonary_lineage.png', dpi=300, bbox_inches='tight')

plt.show()
In [94]:
#this is for the coronary lineage alone exclusively 
score = np.array(terminal_df['Lineage 1'])
y = np.array(terminal_df['C_identity'])

# false positive rate
fpr = []
# true positive rate
tpr = []
# Iterate thresholds from 0.0, 0.01, ... 1.0
thresholds = np.arange(0.0, 1.01, .01)

# get number of positive and negative examples in the dataset
P = sum(y)
N = len(y) - P

# iterate through all thresholds and determine fraction of true positives
# and false positives found at this threshold
for thresh in thresholds:
    FP=0
    TP=0
    for i in range(len(score)):
        if (score[i] > thresh):
            if y[i] == 1:
                TP = TP + 1
            if y[i] == 0:
                FP = FP + 1
    fpr.append(FP/float(N))
    tpr.append(TP/float(P))
    
import sklearn.metrics as metrics
roc_auc = metrics.auc(fpr, tpr)

fig, ax = plt.subplots()
fig.set_figheight(8)
fig.set_figwidth(6)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.title('ROC: Lineage 2 = Coronary Lineage')
plt.plot(fpr, tpr, 'red', label = 'AUC = %0.2f' % roc_auc)
plt.plot([0, 1], [0, 1],'black',linestyle='dashed', label = 'Random Classifier')
plt.legend(loc = 'lower right')
plt.xlim([-0.02, 1])
plt.ylim([0, 1.05])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')

plt.savefig('VSMC/PC_Figures/Control_ROC_coronary_lineage.png', dpi=300, bbox_inches='tight')

plt.show()

Clearly, there is no ability to classify lineage association with tissue identity

This was predicted because Lineage 1 had stationary distance along the Markov transition matrix of 0

Furthermore, less than 1% of cells had an absorption probability associated with lineage 1

Consequently, macrostates and lineages will be re-calculated

In [68]:
g_fwd.compute_macrostates(n_states=1, cluster_key="Tissue")
g_fwd.plot_macrostates(
    discrete=True, legend_loc="right", size=100, basis="fle", title='scCLEAN Macrostates: 1',
    figsize=[8,8], dpi=400)
For 1 macrostate, stationary distribution is computed
Computing eigendecomposition of the transition matrix
Adding `adata.uns['eigendecomposition_fwd']`
       `.eigendecomposition`
    Finish (0:00:02)
Adding `.macrostates`
       `.macrostates_memberships`
       `.coarse_T`
       `.coarse_initial_distribution
       `.coarse_stationary_distribution`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:02)
In [69]:
g_fwd.set_terminal_states_from_macrostates(names=['Pulmonary'])
g_fwd.rename_terminal_states({'Pulmonary':'Lineage_1'})
Adding `adata.obs['terminal_states']`
       `adata.obs['terminal_states_probs']`
       `.terminal_states`
       `.terminal_states_probabilities`
       `.terminal_states_memberships
    Finish`
In [70]:
g_fwd.compute_absorption_probabilities( time_to_absorption='all')
g_fwd.plot_absorption_probabilities(same_plot=False, size=50, basis="fle")
Computing absorption probabilities
WARNING: There is only `1` terminal state, all cells will have probability `1` of going there
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  2.44/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  3.00/s]
Adding `adata.obsm['to_terminal_states']`
       `adata.obsm['absorption_times_fwd']`
       `.absorption_probabilities`
       `.absorption_times`
    Finish (0:00:00)

In [71]:
scv.pl.scatter(adata, basis="fle", c="terminal_states", size=100, legend_loc='on data',
               title='10x-V3 Terminal State', dpi=300, save='VSMC/PC_Figures/Control_1lineage_Terminal_States.png')
saving figure to file VSMC/PC_Figures/Control_1lineage_Terminal_States.png
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: